{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parametric diffusion problem\n",
    "\n",
    " The $\\mu$ is the parameters of the PDE.\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left\\{\\begin{aligned}\n",
    "    -\\nabla \\cdot(\\kappa_\\mu u_\\mu(x)) = f_\\mu(x), & x \\in \\Omega \\\\\n",
    "    u(x) = g_D(x), & x \\in \\Gamma_D \\\\\n",
    "    -\\kappa_\\mu \\nabla u(x) = g_N(x), & x \\in \\Gamma_N \\\\\n",
    "\\end{aligned}\\right.\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_{\\Omega} (u(x) - u(y))\\gamma(x,y) dy   =  f(x)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\int_{\\Omega}u(x)\\gamma(x, y)dy - \\int u(y)\\gamma(x, y)dy = f(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FVM without interpolation.\n",
    "\n",
    "    @article{MIAO2023105,\n",
    "    title = {An interpolation-free cell-centered discretization of the heterogeneous and anisotropic diffusion problems on polygonal meshes},\n",
    "    journal = {Computers & Mathematics with Applications},\n",
    "    volume = {130},\n",
    "    pages = {105-118},\n",
    "    year = {2023},\n",
    "    issn = {0898-1221},\n",
    "    doi = {https://doi.org/10.1016/j.camwa.2022.11.023},\n",
    "    url = {https://www.sciencedirect.com/science/article/pii/S0898122122004928},\n",
    "    author = {Shuai Miao and Jiming Wu and Yanzhong Yao},\n",
    "    keywords = {Interpolation-free, Cell-centered discretization, Coercivity}\n",
    "    }\n",
    "\n",
    " For a squared domain $\\Omega=[a_1, a_2]\\times [b_1, b_2]$. \n",
    " Given a uniform $(N+1)\\times (M+1)$ rectangle mesh $\\mathcal{M} = \\mathcal{v}_{i,j}$ for $i \\in [0, \\dots, N]$ and $j \\in [0, \\dots, M]$. The grid vertexes of mesh $\\mathcal{M}$ should be $\\mathcal{v}_{i, j} = (a_1 + i \\delta_x, ~~b_1 + j \\delta_y)$.\n",
    " One cell $x_{i, j} = ((i-\\frac{1}{2})\\delta_x,~ ~(j-\\frac{1}{2})\\delta_y)$ with size of $\\delta_x = \\frac{(a_2 - a_1)}{N}$, $\\delta_y = \\frac{(b_2 - b_1)}{M}$. \n",
    "\n",
    " Based on the formulas in the paper, we can calculus the construction vector for one cell $x_{i, j}$. \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "    \\nu_{K\\to L} &= x_L - x_K + \\frac{(x_L - x_\\sigma)\\cdot \\mathcal{n}_{K, \\sigma}}{\\lambda^{(n)}_{L, \\sigma}}(\\Lambda_K^T - \\Lambda_L^T)\\mathcal{n}_{K, \\sigma} \\\\\n",
    "    \\\\\n",
    "    \\lambda^{(n)}_{L, \\sigma} &= \\mathcal{n}^T_{K, \\sigma} \\Lambda_K^T \\mathcal{n}_{K, \\sigma}\n",
    "\\end{aligned}\\right.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "\n",
    "<center>\n",
    "\n",
    "<img src=\"./fvm.png\" alt=\"\" width=\"350\" height=\"350\">\n",
    "\n",
    "</center>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the $\\Lambda_K$ is a constant value, then the eq(1) can be reformulated as  the following equation\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\nu_{K\\to L} = x_L - x_K + [(x_L - x_\\sigma)\\cdot n_{K, \\sigma} (\\frac{\\Lambda_K}{\\Lambda_L} - 1)]n_{K, \\sigma}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "And for a uniform rectangle mesh the $(x_L - x_\\sigma)\\cdot n_{K, \\sigma} = \\frac{1}{2}h$ and $x_L - x_K = h n_{K, \\sigma}$. Thus we have\n",
    "\n",
    "$$\n",
    "\\nu_{K\\to L} = h n_{K, \\sigma} + [\\frac{1}{2}h (\\frac{\\Lambda_K}{\\Lambda_L} - 1)]n_{K, \\sigma}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\nu_{K\\to L} = \\frac{3}{2}h (\\frac{\\Lambda_K}{\\Lambda_L} - 1)n_{K, \\sigma}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Given one cel $K = x_{i, j} = ((i-\\frac{1}{2})\\delta_x,~ ~(j-\\frac{1}{2})\\delta_y)$, we can get the following truth:\n",
    " Consider the diffusion matrix is identity matrix. The notation $\\Lambda$ represents the value of diffusion coefficient.\n",
    "<center>\n",
    "\n",
    "$L_{\\kappa}$| $\\sigma_{\\kappa}$| $\\mathcal{n}_{K,\\sigma_{\\kappa}}$ | $x_{L_\\kappa} - \\sigma_\\kappa$ | $x_{L_\\kappa} - x_K$| $\\lambda_{L_\\kappa, \\sigma_\\kappa}^{(n)}$| $\\nu_{K \\to L_{\\kappa}}$|\n",
    "---     | ---    |        ---|        ---|    ---|      ---| ---|\n",
    "$L_1 = ((i-\\frac{3}{2})\\delta_x,~ ~(j-\\frac{1}{2})\\delta_y)^T$ | $\\sigma_1=((i-1)\\delta_x, ~~(j-\\frac{1}{2})\\delta_y)^T$ | $(-1,~~ 0)^T$ | $(-\\frac{1}{2}\\delta_x,~~ 0)^T$| $(-\\delta_x,~~ 0)^T$ | $\\Lambda_K$| $(-\\delta_x,~~ 0)^T + \\frac{1}{2}\\delta_x(1 - \\frac{\\Lambda_{L_1}}{\\Lambda_{K}})(-\\delta_x,~~ 0)^T$ |\n",
    "$L_2 = ((i-\\frac{1}{2})\\delta_x,~ ~(j-\\frac{3}{2})\\delta_y)^T$ | $\\sigma_2=((i-\\frac{1}{2})\\delta_x, ~~(j-1)\\delta_y)^T$ | $(0, ~~-1)^T$ | $(0,~~ -\\frac{1}{2}\\delta_y)^T$| $(0,~~ -\\delta_y)^T$ | $\\Lambda_K$| $(0,~~ -\\delta_y)^T + \\frac{1}{2}\\delta_y(1 - \\frac{\\Lambda_{L_2}}{\\Lambda_{K}})(0,~~ -\\delta_y)^T$ |\n",
    "$L_3 = ((i+\\frac{1}{2})\\delta_x,~ ~(j-\\frac{1}{2})\\delta_y)^T$ | $\\sigma_3=(i\\delta_x, ~~(j-\\frac{1}{2})\\delta_y)^T$     | $(1, ~~0)^T$  |  $(\\frac{1}{2}\\delta_x,~~ 0)^T$|  $(\\delta_x,~~ 0)^T$ | $\\Lambda_K$|  $(\\delta_x,~~ 0)^T + \\frac{1}{2}\\delta_x(1 - \\frac{\\Lambda_{L_3}}{\\Lambda_{K}})(\\delta_x,~~ 0)^T $ |\n",
    "$L_4 = ((i-\\frac{1}{2})\\delta_x,~ ~(j+\\frac{1}{2})\\delta_y)^T$ | $\\sigma_4=((i-\\frac{1}{2})\\delta_x, ~~j\\delta_y)^T$     | $(0, ~~1)^T$  |  $(0,~~ \\frac{1}{2}\\delta_y)^T$|  $(0,~~ \\delta_y)^T$ | $\\Lambda_K$|  $(0,~~ \\delta_y)^T + \\frac{1}{2}\\delta_y(1 - \\frac{\\Lambda_{L_4}}{\\Lambda_{K}})(0,~~ \\delta_y)^T $ |\n",
    "\n",
    "</center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear problems\n",
    "\n",
    " A linear diffusion problem is a PDE, like $\\text{Eq. }(1)$, has a differential operator $\\mathcal{L_{\\mu}} = -\\nabla \\cdot(\\kappa_{\\mu}u_{\\mu})$ which is not related to the unknown function $u$.\n",
    "\n",
    " There are Three problems we mainly want to consider. \n",
    " All of them would be solved on a square domain $\\Omega$ with deep learning methods. \n",
    "\n",
    "### 1. The Heat dissipation of chips\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left\\{\\begin{aligned}\n",
    "    -\\Delta u_\\mu(x) = f_\\mu(x), & x \\in \\Omega \\\\\n",
    "    u(x) = g_D(x), & x \\in \\Gamma_D \\\\\n",
    "    -\\nabla u(x) = g_N(x), & x \\in \\Gamma_N \\\\\n",
    "\\end{aligned}\\right.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "    f(x) = \\sum_{m=1}^M Q_m \\phi_m(x),  \n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "#### Dataset\n",
    "\n",
    " There are $M = 12$ chips in total need to be placed in the domain $\\Omega = [0, 0.1]\\times [0, 0.1]$.\n",
    "\n",
    "<center>\n",
    "\n",
    " ChipsId| Height($m$) | Width($m$)| Intensity($W/m^2$)   |\n",
    " ---    |        ---|        ---|          ---|\n",
    " $\\phi_{1} $ | 0.016     |   0.012   | $Q_{1}  = 4000 $       |  \n",
    " $\\phi_{2} $ | 0.012     |   0.006   | $Q_{2}  = 16000$       |     \n",
    " $\\phi_{3} $ | 0.018     |   0.009   | $Q_{3}  = 6000 $       |  \n",
    " $\\phi_{4} $ | 0.018     |   0.011   | $Q_{4}  = 8000 $       |  \n",
    " $\\phi_{5} $ | 0.018     |   0.018   | $Q_{5}  = 10000$       |     \n",
    " $\\phi_{6} $ | 0.012     |   0.012   | $Q_{6}  = 14000$       |     \n",
    " $\\phi_{7} $ | 0.018     |   0.006   | $Q_{7}  = 16000$       |     \n",
    " $\\phi_{8} $ | 0.009     |   0.009   | $Q_{8}  = 20000$       |     \n",
    " $\\phi_{9} $ | 0.006     |   0.024   | $Q_{9}  = 8000 $       |  \n",
    " $\\phi_{10}$ | 0.006     |   0.012   | $Q_{10} = 16000$       |     \n",
    " $\\phi_{11}$ | 0.012     |   0.024   | $Q_{11} = 10000$       |     \n",
    " $\\phi_{12}$ | 0.024     |   0.024   | $Q_{12} = 20000$       |     \n",
    "\n",
    "</center>\n",
    "\n",
    " Different layouts of Chips can be generated by a random algorithm to construct a dataset which can be regarded as sampling from the random events space. \n",
    " Unlike sampling from a given probability distribution like Gaussian distribution, the distribution of all layouts is unknown in analytic.\n",
    " There are $N_T$ layouts for chips are generated as our dataset for training and $N_V$ layouts for validation. The validation process is used for monitoring the performance of neural networks during the training process. \n",
    "\n",
    "#### Boundary Conditions\n",
    "\n",
    "<center>\n",
    "\n",
    "<img src=\"./chip_bcs.png\" alt=\"\" width=\"750\" height=\"250\">\n",
    "\n",
    "</center>\n",
    "\n",
    " There are three different boundary conditions are considered in out numerical experiments:\n",
    " \n",
    " (a) All sides of $\\Omega$ are Dirichlet boundaries with a const temperature $u = 298K$. \n",
    "\n",
    " (b) The Left side is Dirichlet boundary with const temperature $u = 298K$ and the rest three sides of $\\Omega$ are insulated Neumann boundaries $-\\nabla u = 0$.\n",
    "\n",
    " (c) Only a small sink which is Dirichlet boundary $u = 298K$ on the bottom of $\\Omega$, the rest of $\\partial \\Omega$ are insulated Neumann boundaries $-\\nabla u = 0$.\n",
    " \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.0 The groundwater flow with different point sources\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left\\{\\begin{aligned}\n",
    "    -\\nabla \\cdot \\mathbf{K} u_\\mu(x) = f_\\mu(x), & x \\in \\Omega \\\\\n",
    "    u(x) = g_D(x), & x \\in \\Gamma_D \\\\\n",
    "    -\\nabla u(x) = g_N(x), & x \\in \\Gamma_N \\\\\n",
    "\\end{aligned}\\right.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "    f(x) = \\sum_{m=1}^M Q_m \\delta(x - \\xi_m),  \n",
    "\\end{equation}\n",
    "$$\n",
    " \n",
    " We consider the problem in the domain $\\Omega = [-500, 500]\\times [-500, 500]$. \n",
    " There are ten pump wells at most $M \\in \\{1, 2, \\dots, 10\\}$. \n",
    " The pumping capacity $Q_m$ is limited in the range of $Q_m \\in [5000, 15000] t/d$. \n",
    " The transmissivity tensor $\\mathbf{K} = MK_s$, where $M$ is the thickness of aquifer. We first consider  $ K_s = \\begin{pmatrix} \\frac{100}{3} & 0 \\\\ 0 & \\frac{100}{3} \\end{pmatrix}$ and $M=3$. \n",
    " \n",
    "#### Dataset\n",
    "\n",
    " We first decide how many pumps $M$ randomly. After $M$ is determined, we sample the pumping capacity $Q_m, m=1,\\dots, M$ and its location $\\xi_m$ from the $1d$ uniform distribution in the range of $[5000, 15000]$ and $2d$ uniform distribution in $(-500, 500)\\times (-500, 500)$. \n",
    " There are $N_T$ layouts for pumps are generated as our dataset for training and $N_V$ layouts for validation. The validation process is used for monitoring the performance of neural networks during the training process.\n",
    "\n",
    "#### Boundary Condition\n",
    "\n",
    "<center>\n",
    "\n",
    "<img src=\"./groundwater_bc.png\">\n",
    "\n",
    "</center>\n",
    "\n",
    "The top and bottom sides of $\\Omega$ are Dirichlet boundaries with constant hydraulic head $u = 100$, the left and right sides are non-flux boundaries $-\\nabla u = 0$. \n",
    "\n",
    "### 2.1 The groundwater flow with different point sources\n",
    "\n",
    "We would consider the same equation like above equation but with heterogeneous aquifer, which has a discontinuity at the line $x_1 = 0$. \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "K_s^L = \\begin{pmatrix}\n",
    "    10 & 5 \\\\\n",
    "    5 & 30\n",
    "\\end{pmatrix}, \n",
    "~\n",
    "K_s^R = \\begin{pmatrix}\n",
    "    60 & 10 \\\\\n",
    "    10 & 20\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "<center>\n",
    "\n",
    "<img src=\"./groundwater_bc1.png\">\n",
    "\n",
    "</center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Block-wise diffusion coefficients\n",
    "\n",
    "We consider the following equation with piecewise constant diffusion coefficients. The parameters $\\mu = (\\mu_1, \\dots, \\mu_p) \\in \\mathcal{P} := [0.1,~ 10]^p$. The square domain $\\Omega = [0, 1]\\times [0, 1]$ are decomposed into $p \\in \\{2\\times 2=4,~ 3\\times 3=9,~ 4 \\times 4=16\\}$ pieces $\\Omega_i~, i=1, \\dots, p$. \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left\\{\\begin{aligned}\n",
    "    -\\nabla \\cdot \\kappa_\\mu u_\\mu(x) = 0, &~~ x \\in \\Omega, \\\\\n",
    "    u_\\mu(x) = 0, &~~ x \\in \\Gamma_D, \\\\\n",
    "    -\\nabla u_\\mu(x) = i, &~~ x \\in \\Gamma_{N, i} ~~\\text{for}~~ i = 0, 1,\\\\\n",
    "\\end{aligned}\\right.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\kappa_\\mu(x_1, x_2) := \\sum^p_{i=1} \\mu_i \\mathbb{I}_{\\Omega_i}(x_1, x_2)\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "#### Dataset\n",
    "\n",
    " Different numbers of $M$ are treated as different parametric PDEs. We chose $M=4, 9, 16$ and solve three corresponding parametric PDEs.\n",
    " The dataset was generated by sampled from a uniform distribution with dimension equals to $M$. As we employ the FVM to train our neural networks, There are $N_T + N_V$ matrix need to be assembled.\n",
    "\n",
    "#### Boundary Condition\n",
    "\n",
    " Only the top of the domain $\\Omega$ is Dirichlet boundary with $u_\\mu(x) = 0$ denoted as $\\Gamma_D$. The left and right sides of $\\Omega$ are zero-flux Neumann boundary denoted as $\\Gamma_{N, 0}$. The bottom side is unit flux Neumann boundaru denoted by $\\Gamma_{N, 1}$. \n",
    "\n",
    "<center>\n",
    "<img src=\"./varyk_bc.png\" width=\"300\" height=\"300\">\n",
    "<center>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear problems\n",
    "\n",
    "### 1. \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left\\{\n",
    "  \\begin{aligned}\n",
    "  &-\\nabla\\cdot((u_\\mu + \\mu)\\nabla u) = f, &x \\in \\Omega, \\\\\n",
    "  &u = 0, &x \\in \\partial\\Omega\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "### 2.\n",
    "\n",
    "$$\n",
    "\\kappa(u) = \\left(\\frac{u(1-u)}{u^3 + (1-u)^3}\\right)^2  + \\mu\n",
    "$$\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\left\\{\n",
    "  \\begin{aligned}\n",
    "  &-\\nabla\\cdot(\\kappa(u)\\nabla u) = f, &x \\in \\Omega, \\\\\n",
    "  &u = 0, &x \\in \\partial\\Omega\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For epoch in "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
